#include <fstream>
#include <iostream>
#include <sstream>

#ifndef __CINT__
#include "RooGlobalFunc.h"
#endif
#include "TROOT.h"
#include "TStyle.h"
#include "RooRealVar.h"
#include "RooPolyVar.h"
#include "RooDataSet.h"
#include "RooDataHist.h"
#include "RooAbsData.h"
#include "RooGaussian.h"
#include "RooCBShape.h"
#include "RooAddPdf.h"
#include "RooBifurGauss.h"
#include "RooBreitWigner.h"
#include "RooConstVar.h"
#include "RooErrorVar.h"
#include "RooEffProd.h"
#include "RooChebychev.h"
#include "RooArgusBG.h"
#include "RooPolynomial.h"
#include "RooKeysPdf.h"
#include "RooNDKeysPdf.h"
#include "RooProdPdf.h"
#include "RooFormulaVar.h"
#include "RooHistPdf.h"
#include "RooFitResult.h"
#include "TPaveText.h"
#include "RooRealVar.h"
#include "TCanvas.h"
#include "TLeaf.h"
#include "RooMCStudy.h"
#include "RooRandom.h"
#include "TFile.h"
#include "TH2.h"
#include "TPDF.h"
#include "TRandom3.h"
#include "TChain.h"
#include "TAxis.h"
#include "TH1.h"
#include "RooPlot.h"
#include "RooWorkspace.h"
#include "TFile.h"
#include "RooAddition.h"
#include "RooCategory.h"
#include "RooSimultaneous.h"

using namespace RooFit ;

//
//****************** For Relativistic BW *******************

#include "Riostream.h" 
#include <complex>
#include <vector>
#include "RooRealProxy.h"
#include "RooCategoryProxy.h"
#include "RooAbsCategory.h" 
#include <math.h> 
#include <cmath>
#include "TMath.h" 

class RooRelativisticBreitWigner : public RooAbsPdf {
public:
        RooRelativisticBreitWigner() : m1(0), m2(0) { };
	RooRelativisticBreitWigner(const char *name, const char *title,
			RooAbsReal& _m,
			RooAbsReal& _mass_J,
			RooAbsReal& _width_J,
			RooAbsReal& _r,
			RooAbsReal& _J,
			double _m1,
			double _m2);
	RooRelativisticBreitWigner(const RooRelativisticBreitWigner& other, const char* name=0);
	virtual TObject* clone(const char* newname) const {
		return new RooRelativisticBreitWigner(*this,newname);
	}
	inline virtual ~RooRelativisticBreitWigner() { };
	double CalculateMomentum(const double m) const;
protected:
	RooRealProxy m;
	RooRealProxy mass_J;
	RooRealProxy width_J;
	RooRealProxy r;
	RooRealProxy J;
	const double m1;
	const double m2;
	Double_t evaluate() const;
private:
	ClassDef(RooRelativisticBreitWigner,1) // Your description goes here...
};
//
namespace ComplexPDF {
  double CalculateMomentum(const double m, const double m1, const double m2);
  double BlattWeisskopfFormFactor(const double q, const double q_J, const double r, const unsigned int J);
  std::complex<double> RelativisticBreitWigner(const double m, const double mass_J, const double width_J, const double r, const unsigned int J, const double q, const double q_J);
} // ComplexPDF
//

//**********************************************************
//



//
//************************ utils ***************************

void   util01  (float, int, int, int, int, int, float*, int);
string ftos    (float);
string itos    (int);
float  rounding(float, int );
long   getFSize(string);

//**********************************************************
//



void toy_siml_submit(float ifl, int vnsig, int samset, int order, int bin){

  if(bin<0) return;
  //
  // ROOT Style
  //
  gROOT->Reset();
  gROOT->SetStyle("Plain");
  gStyle->SetOptStat("emr");
  //  gStyle->SetOptStat(1) ;
  gStyle->SetOptTitle(0);
  gStyle->SetPadTopMargin(0.08);
  gStyle->SetPadLeftMargin(0.17);
  gStyle->SetPadBottomMargin(0.13);
  gStyle->SetStatH(0.2);
  gStyle->SetStatW(0.2);
  gStyle->SetTitleSize(0.05,"X");
  gStyle->SetTitleSize(0.05,"Y");
  gStyle->SetTitleOffset(1.00,"X");
  gStyle->SetTitleOffset(1.75,"Y");

  //
  // Options
  //
  int gen    =   500;
  int vnqq   = 19440;
  int vnbb   =   256;
  int vnrare =    70;

  int initNum   = order*bin;
  int lastNum   = (order+1)*bin;
  if(lastNum>gen) lastNum = gen;
  int restNum   = gen%bin;
  int lastOrder = (gen-restNum)/bin;

  //=================== selection number =====================
  //
  //  gsim sig     1
  //  gsim qq     10
  //  gsim bb    100
  //  gsim rare 1000
  //
  //  toy  sig     2
  //  toy  qq     20
  //  toy  bb    200
  //  toy  rare 2000  
  //
  //  --- for example ---
  //
  //  gsim sig  + toy qq   =   21
  //  gsim sig  + toy bb   =  201
  //  gsim sig  + toy rare = 2001
  //  gsim qq   + toy sig  =   12 
  //  gsim qq   + toy bb   =  210 
  //  gsim qq   + toy rare = 2010 
  //  gsim bb   + toy sig  =  102 
  //  gsim bb   + toy qq   =  120 
  //  gsim bb   + toy rare = 2100
  //  gsim rare + toy sig  = 1002  
  //  gsim rare + toy qq   = 1020  
  //  gsim rare + toy bb   = 1200
  //  gsim sig  + toy qq   + gsim bb =  121
  //  gsim sig  + toy qq   + toy  bb =  221
  //
  //==========================================================

  //
  // Save results
  //
  float*   result = new float[60];
  string   output_rt = "result/result_fL"+ftos(ifl)+"_Nsig"+ftos(vnsig)+"_"+itos(samset)+".root";
  ofstream fout;

  TFile file ((output_rt+"_"+itos(order)).c_str(),"recreate");
  TTree *tree = new TTree("h1","h1");

  tree->Branch("v_nsig",   &result[0], "v_nsig/F");
  tree->Branch("e_nsig",   &result[1], "e_nsig/F");
  tree->Branch("v_fl",     &result[2], "v_fl/F");
  tree->Branch("e_fl",     &result[3], "e_fl/F");
  tree->Branch("v_nqq",    &result[4], "v_nqq/F");
  tree->Branch("e_nqq",    &result[5], "e_nqq/F");
  tree->Branch("v_nbb",    &result[6], "v_nbb/F");
  tree->Branch("e_nbb",    &result[7], "e_nbb/F");
  tree->Branch("v_deqqa0", &result[8], "v_deqqa0/F");
  tree->Branch("e_deqqa0", &result[9], "e_deqqa0/F");
  tree->Branch("v_mbcqqp", &result[10], "v_mbcqqp/F");
  tree->Branch("e_mbcqqp", &result[11], "e_mbcqqp/F");
  tree->Branch("v_k0mqq0", &result[12], "v_k0mqq0/F");
  tree->Branch("e_k0mqq0", &result[13], "e_k0mqq0/F");
  tree->Branch("v_k0mqqm", &result[14], "v_k0mqqm/F");
  tree->Branch("e_k0mqqm", &result[15], "e_k0mqqm/F");
  tree->Branch("v_k0mqqw", &result[16], "v_k0mqqw/F");
  tree->Branch("e_k0mqqw", &result[17], "e_k0mqqw/F");
  tree->Branch("v_k0mqqf", &result[18], "v_k0mqqf/F");
  tree->Branch("e_k0mqqf", &result[19], "e_k0mqqf/F");
  tree->Branch("v_kpmqq0", &result[20], "v_kpmqq0/F");
  tree->Branch("e_kpmqq0", &result[21], "e_kpmqq0/F");
  tree->Branch("v_kpmqqm", &result[22], "v_kpmqqm/F");
  tree->Branch("e_kpmqqm", &result[23], "e_kpmqqm/F");
  tree->Branch("v_kpmqqw", &result[24], "v_kpmqqw/F");
  tree->Branch("e_kpmqqw", &result[25], "e_kpmqqw/F");
  tree->Branch("v_kpmqqf", &result[26], "v_kpmqqf/F");
  tree->Branch("e_kpmqqf", &result[27], "e_kpmqqf/F");
  tree->Branch("v_nnqqgm", &result[28], "v_nnqqgm/F");
  tree->Branch("e_nnqqgm", &result[29], "e_nnqqgm/F");
  tree->Branch("v_nnqqgs", &result[30], "v_nnqqgs/F");
  tree->Branch("e_nnqqgs", &result[31], "e_nnqqgs/F");

  tree->Branch("v_k0hcs0", &result[32], "v_k0hcs0/F");
  tree->Branch("v_k0hcs1", &result[33], "v_k0hcs1/F");
  tree->Branch("v_k0hcs2", &result[34], "v_k0hcs2/F");
  tree->Branch("v_k0hcs3", &result[35], "v_k0hcs3/F");

  tree->Branch("v_kphcs0", &result[36], "v_kphcs0/F");
  tree->Branch("v_kphcs1", &result[37], "v_kphcs1/F");
  tree->Branch("v_kphcs2", &result[38], "v_kphcs2/F");
  tree->Branch("v_kphcs3", &result[39], "v_kphcs3/F");
  tree->Branch("v_kphcs4", &result[40], "v_kphcs4/F");

  tree->Branch("e_k0hcs0", &result[41], "e_k0hcs0/F");
  tree->Branch("e_k0hcs1", &result[42], "e_k0hcs1/F");
  tree->Branch("e_k0hcs2", &result[43], "e_k0hcs2/F");
  tree->Branch("e_k0hcs3", &result[44], "e_k0hcs3/F");

  tree->Branch("e_kphcs0", &result[45], "e_kphcs0/F");
  tree->Branch("e_kphcs1", &result[46], "e_kphcs1/F");
  tree->Branch("e_kphcs2", &result[47], "e_kphcs2/F");
  tree->Branch("e_kphcs3", &result[48], "e_kphcs3/F");
  tree->Branch("e_kphcs4", &result[49], "e_kphcs4/F");

  tree->Branch("v_debb0",  &result[50], "v_debb0/F");
  tree->Branch("v_mbcbb0", &result[51], "v_mbcbb0/F");
  tree->Branch("v_k0mbb0", &result[52], "v_k0mbb0/F");
  tree->Branch("v_kpmbb0", &result[53], "v_kpmbb0/F");

  tree->Branch("v_nnqqgf", &result[54], "v_nnqqgf/F");
  tree->Branch("e_nnqqgf", &result[55], "e_nnqqgf/F");

  for(int i=initNum; i<lastNum; i++){

    cout<<"*************************************************************************************************"<<endl;
    cout<<""<<endl;    cout<<""<<endl;    cout<<""<<endl;    cout<<""<<endl;    cout<<""<<endl;
    cout<<"                                fL = "<<ifl<<" & Nsig = "<<vnsig<<"  ->  "<<i+1<<"               "<<endl;
    cout<<""<<endl;    cout<<""<<endl;    cout<<""<<endl;    cout<<""<<endl;    cout<<""<<endl;
    cout<<"*************************************************************************************************"<<endl;

    util01( ifl, vnsig, vnqq, vnbb, vnrare, i, result, samset );

    tree->Fill();

  }

  tree->GetCurrentFile()->Write();
  file.Close();

  string addOutput = "hadd " + output_rt;     
  string rmOutput  = "rm -f ";

  for(int i=0; i<lastOrder; i++){ 
    addOutput = addOutput + " "  + output_rt + "_" + itos(i);         
    rmOutput  = rmOutput  + " "  + output_rt + "_" + itos(i);         
  }

  // check filesize
  float fileSize = 0.0;
  float tmpSize  = 0.0;
  for(int i=0; i<lastOrder; i++){
    tmpSize = getFSize((output_rt + "_" + itos(i)).c_str());
    if(fileSize<tmpSize){
      fileSize = tmpSize;
    }
  }
  cout<<fileSize<<endl;
          
  int checkSize = 1;  
  for(int i=0; i<lastOrder; i++){ 
    if(getFSize((output_rt + "_" + itos(i)).c_str())<fileSize*0.99) checkSize = 0;
  }
          
  if(checkSize>0){   
    system(("rm -f "+output_rt).c_str());     
    system(addOutput.c_str());    
    system(rmOutput.c_str());     
  }

  return;
}




void util01(float ifl, int vnsig, int vnqq, int vnbb, int vnrare, int num, float* result, int samset){

  //
  // Call each files for each experiment
  //
  string in_root1  = "../mk_sample/signals/sig_fL"+ ftos(ifl) +"_Nsig" + itos(vnsig) + "_" + itos(num) +".root";
  if(ifl==0.75){
    cout<<"fL=0.75 signal MC samples are included !!! "<<endl;
    in_root1  = "../mk_sample/signal_fL075/sig_fL"+ ftos(ifl) +"_Nsig" + itos(vnsig) + "_" + itos(num) +".root";
  }
  string in_root2  = "../mk_sample/samples/qq_"   + itos(num) +".root";
  string in_root3  = "../mk_sample/samples/bb_"   + itos(num) +".root";
  string in_root4  = "../mk_sample/samples/rare_" + itos(num) +".root";

  //
  // Set numbers for each MC type 
  //   
  //----------------------------------------
  //  1: gsim   2:toy    0:none
  // 
  //  gsim : sig=1  qq=10  bb=100 rare=1000
  //  toy  : sig=2  qq=20  bb=200 rare=2000
  //----------------------------------------
  int samset01 = (int) samset;
  int samset02 = (int) samset/10;
  int samset03 = (int) samset/100;
  int samset04 = (int) samset/1000;

  int set_sig  = samset01%5; 
  int set_qq   = samset02%5; 
  int set_bb   = samset03%5;
  int set_rare = samset04%5;

  //
  //********************************** fit start ************************************

  //
  // true signal fraction
  //
  Float_t true_frac1 = 1.0 - 0.07742;
  Float_t true_frac0 = 1.0 - 0.03518;
  Float_t true_frac1A = 1.0 - 0.2015;
  Float_t true_frac1B = 1.0 - 0.0716;
  Float_t true_frac0A = 1.0 - 0.0709;
  Float_t true_frac0B = 1.0 - 0.0333;

  RooRealVar true_fraction1("true_fraction1", "lo true signal fraction", true_frac1);
  RooRealVar true_fraction0("true_fraction0", "tr true signal fraction", true_frac0);
  RooRealVar true_fraction1A("true_fraction1A", "lo true signal fraction", true_frac1A);
  RooRealVar true_fraction1B("true_fraction1B", "lo true signal fraction", true_frac1B);
  RooRealVar true_fraction0A("true_fraction0A", "tr true signal fraction", true_frac0A);
  RooRealVar true_fraction0B("true_fraction0B", "tr true signal fraction", true_frac0B);                                                                  


  //
  // Define Variables
  //
  RooRealVar de      ("de",       "#DeltaE",     -0.15, 0.15  );
  RooRealVar mbc     ("mbc",       "M_{bc}",      5.25, 5.2889);
  RooRealVar k0_mass ("k0_mass",   "mass(K*0)",   0.78, 1.00  );
  RooRealVar kp_mass ("kp_mass",   "mass(K*+)",   0.78, 1.00  );
  RooRealVar k0_hcos ("k0_hcos", "K*0 helicity", -1.00, 1.00  );
  RooRealVar kp_hcos ("kp_hcos", "K*+ helicity", -1.00, 1.00  );
  RooRealVar nn_tran2("nn_tran2", "N_{trans}",  -12.0, 12.0   );

  //
  // Variable Set
  //
  RooArgList *varset = new RooArgList(de, mbc, k0_mass, kp_mass, k0_hcos, kp_hcos, nn_tran2);

  //
  // RooCategory
  //
  RooCategory sample("sample","sample") ;
  sample.defineType("r01") ;
  sample.defineType("r02") ;

  //
  // Call HistPdf info
  //
  TFile *f11 = new TFile("histpdf/de_scf_lo.root");
  TFile *f12 = new TFile("histpdf/de_scf_tr.root");
  TFile *f21 = new TFile("histpdf/mbc_scf_lo.root");
  TFile *f22 = new TFile("histpdf/mbc_scf_tr.root");
  TFile *f31 = new TFile("histpdf/k0mass_scf_lo.root");
  TFile *f32 = new TFile("histpdf/k0mass_scf_tr.root");
  TFile *f41 = new TFile("histpdf/kpmass_scf_lo.root");
  TFile *f42 = new TFile("histpdf/kpmass_scf_tr.root");

  TFile *f51 = new TFile("histpdf/k0hcos_true_lo01.root");
  TFile *f52 = new TFile("histpdf/k0hcos_true_lo02.root");
  TFile *f53 = new TFile("histpdf/k0hcos_true_tr01.root");
  TFile *f54 = new TFile("histpdf/k0hcos_true_tr02.root");
  TFile *f55 = new TFile("histpdf/k0hcos_scf_lo01.root");
  TFile *f56 = new TFile("histpdf/k0hcos_scf_lo02.root");
  TFile *f57 = new TFile("histpdf/k0hcos_scf_tr01.root");
  TFile *f58 = new TFile("histpdf/k0hcos_scf_tr02.root");
  TFile *f59 = new TFile("histpdf/k0hcos_rare01.root");
  TFile *f60 = new TFile("histpdf/k0hcos_rare02.root");

  TFile *f61 = new TFile("histpdf/kphcos_true_lo.root");
  TFile *f62 = new TFile("histpdf/kphcos_true_tr.root");
  TFile *f63 = new TFile("histpdf/kphcos_scf_lo.root");
  TFile *f64 = new TFile("histpdf/kphcos_scf_tr.root");
  TFile *f67 = new TFile("histpdf/kphcos_rare.root");

  RooWorkspace* w11 = (RooWorkspace*) f11->Get("w");
  RooWorkspace* w12 = (RooWorkspace*) f12->Get("w");
  RooWorkspace* w21 = (RooWorkspace*) f21->Get("w");
  RooWorkspace* w22 = (RooWorkspace*) f22->Get("w");
  RooWorkspace* w31 = (RooWorkspace*) f31->Get("w");
  RooWorkspace* w32 = (RooWorkspace*) f32->Get("w");
  RooWorkspace* w41 = (RooWorkspace*) f41->Get("w");
  RooWorkspace* w42 = (RooWorkspace*) f42->Get("w");
  RooWorkspace* w51 = (RooWorkspace*) f51->Get("w");
  RooWorkspace* w52 = (RooWorkspace*) f52->Get("w");
  RooWorkspace* w53 = (RooWorkspace*) f53->Get("w");
  RooWorkspace* w54 = (RooWorkspace*) f54->Get("w");
  RooWorkspace* w55 = (RooWorkspace*) f55->Get("w");
  RooWorkspace* w56 = (RooWorkspace*) f56->Get("w");
  RooWorkspace* w57 = (RooWorkspace*) f57->Get("w");
  RooWorkspace* w58 = (RooWorkspace*) f58->Get("w");
  RooWorkspace* w59 = (RooWorkspace*) f59->Get("w");
  RooWorkspace* w60 = (RooWorkspace*) f60->Get("w");
  RooWorkspace* w61 = (RooWorkspace*) f61->Get("w");
  RooWorkspace* w62 = (RooWorkspace*) f62->Get("w");
  RooWorkspace* w63 = (RooWorkspace*) f63->Get("w");
  RooWorkspace* w64 = (RooWorkspace*) f64->Get("w");
  RooWorkspace* w67 = (RooWorkspace*) f67->Get("w");

  RooAbsPdf* model11 = w11->pdf("de_scf_lo_pdf");
  RooAbsPdf* model12 = w12->pdf("de_scf_tr_pdf");
  RooAbsPdf* model21 = w21->pdf("mbc_scf_lo_pdf");
  RooAbsPdf* model22 = w22->pdf("mbc_scf_tr_pdf");
  RooAbsPdf* model31 = w31->pdf("k0mass_scf_lo_pdf");
  RooAbsPdf* model32 = w32->pdf("k0mass_scf_tr_pdf");
  RooAbsPdf* model41 = w41->pdf("kpmass_scf_lo_pdf");
  RooAbsPdf* model42 = w42->pdf("kpmass_scf_tr_pdf");

  RooAbsPdf* model51 = w51->pdf("k0hcos_true_lo_pdf01");
  RooAbsPdf* model52 = w52->pdf("k0hcos_true_lo_pdf02");
  RooAbsPdf* model53 = w53->pdf("k0hcos_true_tr_pdf01");
  RooAbsPdf* model54 = w54->pdf("k0hcos_true_tr_pdf02");
  RooAbsPdf* model55 = w55->pdf("k0hcos_scf_lo_pdf01");
  RooAbsPdf* model56 = w56->pdf("k0hcos_scf_lo_pdf02");
  RooAbsPdf* model57 = w57->pdf("k0hcos_scf_tr_pdf01");
  RooAbsPdf* model58 = w58->pdf("k0hcos_scf_tr_pdf02");
  RooAbsPdf* model59 = w59->pdf("k0hcos_rare_pdf01");
  RooAbsPdf* model60 = w60->pdf("k0hcos_rare_pdf02");

  RooAbsPdf* model61 = w61->pdf("kphcos_true_lo_pdf");
  RooAbsPdf* model62 = w62->pdf("kphcos_true_tr_pdf");
  RooAbsPdf* model63 = w63->pdf("kphcos_scf_lo_pdf");
  RooAbsPdf* model64 = w64->pdf("kphcos_scf_tr_pdf");
  RooAbsPdf* model67 = w67->pdf("kphcos_rare_pdf");

  //
  // Signal PDF Setting
  //

  // parameters for relativistic BW
  const double m_ks = 0.497614;
  const double m_pi = 0.13957018;
  const double m_ka = 0.493677;
  RooRealVar spin ("spin", "spin", 1.0);

  //-------------------------------------- de -----------------------------------------------

  //
  // true de
  //
  /*
  // Gaussian1
  RooRealVar de_true_sig_g1_mean0("de_true_sig_g1_mean0", "de true signal mean",   8.4e-04 );
  RooRealVar de_true_sig_g1_mean1("de_true_sig_g1_mean1", "de true signal mean",   7.0e-04 );
  RooRealVar de_true_sig_g1_wid0 ("de_true_sig_g1_wid0",  "de true signal width",  0.00996 );
  RooRealVar de_true_sig_g1_wid1 ("de_true_sig_g1_wid1",  "de true signal width",  0.01086 );

  // Gaussian2
  RooRealVar de_true_sig_g2_mean0("de_true_sig_g2_mean0", "de true signal mean", -0.00868 );
  RooRealVar de_true_sig_g2_mean1("de_true_sig_g2_mean1", "de true signal mean", -0.00900 );
  RooRealVar de_true_sig_g2_wid0 ("de_true_sig_g2_wid0",  "de true signal width", 0.04861 );
  RooRealVar de_true_sig_g2_wid1 ("de_true_sig_g2_wid1",  "de true signal width", 0.04909 );

  // fraction of Gaussian1
  RooRealVar de_true_sig_frac0("de_true_sig_frac0", "de true signal fraction", 0.9017 );
  RooRealVar de_true_sig_frac1("de_true_sig_frac1", "de true signal fraction", 0.8987 );
  */
  // Gaussian1
  RooRealVar de_true_sig_g1_mean0("de_true_sig_g1_mean0", "de true signal mean",   8.0e-04 );
  RooRealVar de_true_sig_g1_mean1("de_true_sig_g1_mean1", "de true signal mean",   6.6e-04 );
  RooRealVar de_true_sig_g1_wid0 ("de_true_sig_g1_wid0",  "de true signal width",  0.00896 );
  RooRealVar de_true_sig_g1_wid1 ("de_true_sig_g1_wid1",  "de true signal width",  0.00955 );

  // Gaussian2
  RooRealVar de_true_sig_g2_mean0("de_true_sig_g2_mean0", "de true signal mean",  0.00160 );
  RooRealVar de_true_sig_g2_mean1("de_true_sig_g2_mean1", "de true signal mean", -0.00004 );
  RooRealVar de_true_sig_g2_wid0 ("de_true_sig_g2_wid0",  "de true signal width", 0.00896 );
  RooRealVar de_true_sig_g2_wid1 ("de_true_sig_g2_wid1",  "de true signal width", 0.02241 );

  // fraction of Gaussian1
  RooRealVar de_true_sig_frac0("de_true_sig_frac0", "de true signal fraction", 0.8026 );
  RooRealVar de_true_sig_frac1("de_true_sig_frac1", "de true signal fraction", 0.7647 );

  RooGaussian de_true_sig_g1_pdf0("de_true_sig_g1_pdf0", "de true signal Gaussian", de, de_true_sig_g1_mean0, de_true_sig_g1_wid0 );
  RooGaussian de_true_sig_g1_pdf1("de_true_sig_g1_pdf1", "de true signal Gaussian", de, de_true_sig_g1_mean1, de_true_sig_g1_wid1 );
  RooGaussian de_true_sig_g2_pdf0("de_true_sig_g2_pdf0", "de true signal Gaussian", de, de_true_sig_g2_mean0, de_true_sig_g2_wid0 );
  RooGaussian de_true_sig_g2_pdf1("de_true_sig_g2_pdf1", "de true signal Gaussian", de, de_true_sig_g2_mean1, de_true_sig_g2_wid1 );

  // Gaussian1 + Gaussian2
  RooAddPdf de_true_sig_pdf0("de_true_sig_pdf0", "de true signal pdf", RooArgList(de_true_sig_g1_pdf0, de_true_sig_g2_pdf0), de_true_sig_frac0);
  RooAddPdf de_true_sig_pdf1("de_true_sig_pdf1", "de true signal pdf", RooArgList(de_true_sig_g1_pdf1, de_true_sig_g2_pdf1), de_true_sig_frac1);

  //----------------------------------- mbc ----------------------------------------------

  // Gaussian 
  RooRealVar mbc_true_sig_gau_mean0("mbc_true_sig_gau_mean0", "mbc true signal mean",  5.27905);
  RooRealVar mbc_true_sig_gau_mean1("mbc_true_sig_gau_mean1", "mbc true signal mean",  5.27904);
  RooRealVar mbc_true_sig_gau_wid0 ("mbc_true_sig_gau_wid0",  "mbc true signal width", 0.002609 );
  RooRealVar mbc_true_sig_gau_wid1 ("mbc_true_sig_gau_wid1",  "mbc true signal width", 0.002615 );

  RooGaussian mbc_true_sig_pdf0("mbc_true_sig_pdf0",  "mbc true signal Gaussian", mbc, mbc_true_sig_gau_mean0, mbc_true_sig_gau_wid0);
  RooGaussian mbc_true_sig_pdf1("mbc_true_sig_pdf1",  "mbc true signal Gaussian", mbc, mbc_true_sig_gau_mean1, mbc_true_sig_gau_wid1);

  //----------------------------------- k0 mass ---------------------------------------------

  // relativistic BW 
  RooRealVar k0mass_true_sig_bw_mean0("k0mass_true_sig_bw_mean0", "k0mass true signal mean",  0.895109);
  RooRealVar k0mass_true_sig_bw_mean1("k0mass_true_sig_bw_mean1", "k0mass true signal mean",  0.894144);
  RooRealVar k0mass_true_sig_bw_wid0 ("k0mass_true_sig_bw_wid0",  "k0mass true signal width", 0.05160);
  RooRealVar k0mass_true_sig_bw_wid1 ("k0mass_true_sig_bw_wid1",  "k0mass true signal width", 0.05025);
  RooRealVar k0mass_true_sig_bw_par0 ("k0mass_true_sig_bw_par0",  "k0mass true signal par",   5.557);
  RooRealVar k0mass_true_sig_bw_par1 ("k0mass_true_sig_bw_par1",  "k0mass true signal par",   4.116);

  RooRelativisticBreitWigner k0mass_true_sig_pdf0("k0mass_true_sig_pdf0",  "k0mass true signal BW", k0_mass, k0mass_true_sig_bw_mean0, k0mass_true_sig_bw_wid0, k0mass_true_sig_bw_par0, spin, m_pi, m_ka);
  RooRelativisticBreitWigner k0mass_true_sig_pdf1("k0mass_true_sig_pdf1",  "k0mass true signal BW", k0_mass, k0mass_true_sig_bw_mean1, k0mass_true_sig_bw_wid1, k0mass_true_sig_bw_par1, spin, m_pi, m_ka);

  //----------------------------------- kp mass ---------------------------------------------

  // relativistic BW 
  RooRealVar kpmass_true_sig_bw_mean0("kpmass_true_sig_bw_mean0", "kpmass true signal mean",  0.890777);
  RooRealVar kpmass_true_sig_bw_mean1("kpmass_true_sig_bw_mean1", "kpmass true signal mean",  0.890635);
  RooRealVar kpmass_true_sig_bw_wid0 ("kpmass_true_sig_bw_wid0",  "kpmass true signal width", 0.05207);
  RooRealVar kpmass_true_sig_bw_wid1 ("kpmass_true_sig_bw_wid1",  "kpmass true signal width", 0.05218);
  RooRealVar kpmass_true_sig_bw_par0 ("kpmass_true_sig_bw_par0",  "kpmass true signal par",   5.850);
  RooRealVar kpmass_true_sig_bw_par1 ("kpmass_true_sig_bw_par1",  "kpmass true signal par",   6.460);

  RooRelativisticBreitWigner kpmass_true_sig_pdf0("kpmass_true_sig_pdf0",  "kpmass true signal BW", kp_mass, kpmass_true_sig_bw_mean0, kpmass_true_sig_bw_wid0, kpmass_true_sig_bw_par0, spin, m_ks, m_pi);
  RooRelativisticBreitWigner kpmass_true_sig_pdf1("kpmass_true_sig_pdf1",  "kpmass true signal BW", kp_mass, kpmass_true_sig_bw_mean1, kpmass_true_sig_bw_wid1, kpmass_true_sig_bw_par1, spin, m_ks, m_pi);

  //----------------------------------- k0 hcos ---------------------------------------------

                     /************* HistPdf modelling **************/

  RooSimultaneous k0hcos_true_sig_pdf0("k0hcos_true_sig_pdf0","simultaneous pdf",sample) ;
  k0hcos_true_sig_pdf0.addPdf(*model53, "r01");
  k0hcos_true_sig_pdf0.addPdf(*model54, "r02");

  RooSimultaneous k0hcos_true_sig_pdf1("k0hcos_true_sig_pdf1","simultaneous pdf",sample) ;
  k0hcos_true_sig_pdf1.addPdf(*model51, "r01");
  k0hcos_true_sig_pdf1.addPdf(*model52, "r02");

  RooSimultaneous k0hcos_scf_sig_pdf0("k0hcos_scf_sig_pdf0","simultaneous pdf",sample) ;
  k0hcos_scf_sig_pdf0.addPdf(*model57, "r01");
  k0hcos_scf_sig_pdf0.addPdf(*model58, "r02");

  RooSimultaneous k0hcos_scf_sig_pdf1("k0hcos_scf_sig_pdf1","simultaneous pdf",sample) ;
  k0hcos_scf_sig_pdf1.addPdf(*model55, "r01");
  k0hcos_scf_sig_pdf1.addPdf(*model56, "r02");

  //----------------------------------- kp hcos ---------------------------------------------

                      /************* HistPdf modelling **************/

  //----------------------------------- NN out ----------------------------------------------

  // BifurGauss1 
  RooRealVar nn_tran2_true_sig_bif1_mean0 ("nn_tran2_true_sig_bif1_mean0", "nn_tran2 true signal mean",     4.8100);
  RooRealVar nn_tran2_true_sig_bif1_widL0 ("nn_tran2_true_sig_bif1_widL0",  "nn_tran2 true signal width L", 3.4500);
  RooRealVar nn_tran2_true_sig_bif1_widR0 ("nn_tran2_true_sig_bif1_widR0",  "nn_tran2 true signal width R", 0.8400);
  RooRealVar nn_tran2_true_sig_bif1_mean1 ("nn_tran2_true_sig_bif1_mean1", "nn_tran2 true signal mean",     4.7370);
  RooRealVar nn_tran2_true_sig_bif1_widL1 ("nn_tran2_true_sig_bif1_widL1",  "nn_tran2 true signal width L", 3.4210);
  RooRealVar nn_tran2_true_sig_bif1_widR1 ("nn_tran2_true_sig_bif1_widR1",  "nn_tran2 true signal width R", 0.8570);

  // BifurGauss2 
  RooRealVar nn_tran2_true_sig_bif2_mean0 ("nn_tran2_true_sig_bif2_mean0", "nn_tran2 true signal mean",     3.2880);
  RooRealVar nn_tran2_true_sig_bif2_widL0 ("nn_tran2_true_sig_bif2_widL0",  "nn_tran2 true signal width L", 2.2330);
  RooRealVar nn_tran2_true_sig_bif2_widR0 ("nn_tran2_true_sig_bif2_widR0",  "nn_tran2 true signal width R", 1.0700);
  RooRealVar nn_tran2_true_sig_bif2_mean1 ("nn_tran2_true_sig_bif2_mean1", "nn_tran2 true signal mean",     3.2660);
  RooRealVar nn_tran2_true_sig_bif2_widL1 ("nn_tran2_true_sig_bif2_widL1",  "nn_tran2 true signal width L", 2.1960);
  RooRealVar nn_tran2_true_sig_bif2_widR1 ("nn_tran2_true_sig_bif2_widR1",  "nn_tran2 true signal width R", 1.0790);

  // fraction between BifurGauss1 + BifurGauss2
  RooRealVar nn_tran2_true_sig_frac0("nn_tran2_true_sig_frac0", "nn_tran2 true signal fraction 2 Bifurs", 0.3000);
  RooRealVar nn_tran2_true_sig_frac1("nn_tran2_true_sig_frac1", "nn_tran2 true signal fraction 2 Bifurs", 0.3310);

  // BifurGauss (SCF)
  RooRealVar nn_tran2_scf_sig_bif1_mean0 ("nn_tran2_scf_sig_bif1_mean0",  "nn_tran2 scf signal mean",    2.8720);
  RooRealVar nn_tran2_scf_sig_bif1_widL0 ("nn_tran2_scf_sig_bif1_widL0",  "nn_tran2 scf signal width L", 2.5580);
  RooRealVar nn_tran2_scf_sig_bif1_widR0 ("nn_tran2_scf_sig_bif1_widR0",  "nn_tran2 scf signal width R", 1.3540);
  RooRealVar nn_tran2_scf_sig_bif1_mean1 ("nn_tran2_scf_sig_bif1_mean1",  "nn_tran2 scf signal mean",    3.0100);
  RooRealVar nn_tran2_scf_sig_bif1_widL1 ("nn_tran2_scf_sig_bif1_widL1",  "nn_tran2 scf signal width L", 2.5160);
  RooRealVar nn_tran2_scf_sig_bif1_widR1 ("nn_tran2_scf_sig_bif1_widR1",  "nn_tran2 scf signal width R", 1.2810);

  RooBifurGauss nn_tran2_true_sig_bif1_pdf0("nn_tran2_true_sig_bif1_pdf0",  "nn_tran2 true signal BifurGauss", nn_tran2, nn_tran2_true_sig_bif1_mean0, nn_tran2_true_sig_bif1_widL0, nn_tran2_true_sig_bif1_widR0);
  RooBifurGauss nn_tran2_true_sig_bif1_pdf1("nn_tran2_true_sig_bif1_pdf1",  "nn_tran2 true signal BifurGauss", nn_tran2, nn_tran2_true_sig_bif1_mean1, nn_tran2_true_sig_bif1_widL1, nn_tran2_true_sig_bif1_widR1);
  RooBifurGauss nn_tran2_true_sig_bif2_pdf0("nn_tran2_true_sig_bif2_pdf0",  "nn_tran2 true signal BifurGauss", nn_tran2, nn_tran2_true_sig_bif2_mean0, nn_tran2_true_sig_bif2_widL0, nn_tran2_true_sig_bif2_widR0);
  RooBifurGauss nn_tran2_true_sig_bif2_pdf1("nn_tran2_true_sig_bif2_pdf1",  "nn_tran2 true signal BifurGauss", nn_tran2, nn_tran2_true_sig_bif2_mean1, nn_tran2_true_sig_bif2_widL1, nn_tran2_true_sig_bif2_widR1);

  // BifurGauss1 + BifurGauss2 
  RooAddPdf nn_tran2_true_sig_pdf0("nn_tran2_true_sig_pdf0", "nn_tran2 true signal pdf", RooArgList(nn_tran2_true_sig_bif1_pdf0, nn_tran2_true_sig_bif2_pdf0), nn_tran2_true_sig_frac0);
  RooAddPdf nn_tran2_true_sig_pdf1("nn_tran2_true_sig_pdf1", "nn_tran2 true signal pdf", RooArgList(nn_tran2_true_sig_bif1_pdf1, nn_tran2_true_sig_bif2_pdf1), nn_tran2_true_sig_frac1);

  // BifurGauss (SCF)
  RooBifurGauss nn_tran2_scf_sig_pdf0("nn_tran2_scf_sig_pdf0",  "nn_tran2 scf signal BifurGauss", nn_tran2, nn_tran2_scf_sig_bif1_mean0, nn_tran2_scf_sig_bif1_widL0, nn_tran2_scf_sig_bif1_widR0);
  RooBifurGauss nn_tran2_scf_sig_pdf1("nn_tran2_scf_sig_pdf1",  "nn_tran2 scf signal BifurGauss", nn_tran2, nn_tran2_scf_sig_bif1_mean1, nn_tran2_scf_sig_bif1_widL1, nn_tran2_scf_sig_bif1_widR1);

  //
  // qq PDF Setting
  //

  //-------------------------------------- de -----------------------------------------------

  // Cheby 
  RooRealVar   de_qq_cheby("a0_{Cheby}^{#DeltaE}", "a_{0}(#DeltaE,Cheby)", -0.15328, -0.3, 0.0);

  RooChebychev de_qq_pdf("de_qq_pdf" , "de qq pdf", de, RooArgList(de_qq_cheby));


  //----------------------------------- mbc ----------------------------------------------

  // Argus 
  RooRealVar mbc_qq_argus_cutoff("mbc_qq_argus_cutoff", "cutoff(#M_{bc}^{new},Argus)",   5.2889 ); 
  RooRealVar mbc_qq_argus_par   ("shape_{Argus}",          "shape(#M_{bc}^{new},Argus)",  -30.483, -70.0, -5.0 );

  RooArgusBG mbc_qq_pdf( "mbc_qq_pdf", "mbc qq pdf", mbc, mbc_qq_argus_cutoff, mbc_qq_argus_par);

  //----------------------------------- k0 mass ---------------------------------------------

  // relativistic BW
  RooRealVar k0mass_qq_bw_mean("k0mass_qq_bw_mean", "k0mass qq mean",  0.89575);
  RooRealVar k0mass_qq_bw_wid ("k0mass_qq_bw_wid",  "k0mass qq width", 0.0586);
  RooRealVar k0mass_qq_bw_par ("k0mass_qq_bw_par",  "k0mass qq par",   6.46);

  // Cheby 
  RooRealVar k0mass_qq_cheby01("k0mass_qq_cheby01",  "a_{0}(m(K*0),Cheby)",  0.0328, -0.3, 0.3);

  // fraction between relativistic BW + Cheby
  RooRealVar bw_qq_fraction1  ("bw_qq_fraction1", "bw qq fraction", 0.3461, 0.0, 0.7);

  RooRelativisticBreitWigner k0mass_qq_bw_pdf   ("k0mass_qq_bw_pdf",    "k0mass qq BW",         k0_mass, k0mass_qq_bw_mean, k0mass_qq_bw_wid, k0mass_qq_bw_par, spin, m_pi, m_ka );
  RooChebychev               k0mass_qq_cheby_pdf("k0mass_qq_cheby_pdf", "k0_mass qq cheby pdf", k0_mass, RooArgList(k0mass_qq_cheby01));

  // relativistic BW + Cheby
  RooAddPdf k0mass_qq_pdf("k0mass_qq_pdf", "k0mass qq pdf", RooArgList(k0mass_qq_bw_pdf, k0mass_qq_cheby_pdf), bw_qq_fraction1);


  //----------------------------------- kp mass ---------------------------------------------

  // relativistic BW
  RooRealVar kpmass_qq_bw_mean("kpmass_qq_bw_mean", "kpmass qq mean",  0.89222);
  RooRealVar kpmass_qq_bw_wid ("kpmass_qq_bw_wid",  "kpmass qq width", 0.0554);
  RooRealVar kpmass_qq_bw_par ("kpmass_qq_bw_par",  "kpmass qq par",   8.19);

  // Cheby
  RooRealVar   kpmass_qq_cheby01( "a0_{Cheby}^{K*+}","a_{0}(m(K*+),Cheby)", -0.02699, -0.3, 0.3);

  // fraction between relativistic BW + Cheby
  RooRealVar bw_qq_fraction2("f_{BW}^{K*+}", "bw qq fraction", 0.4022, 0.0, 0.8);

  RooRelativisticBreitWigner kpmass_qq_bw_pdf   ("kpmass_qq_bw_pdf",    "kpmass qq BW",         kp_mass, kpmass_qq_bw_mean, kpmass_qq_bw_wid, kpmass_qq_bw_par, spin, m_ks, m_pi );  
  RooChebychev               kpmass_qq_cheby_pdf("kpmass_qq_cheby_pdf", "kp_mass qq cheby pdf", kp_mass, RooArgList(kpmass_qq_cheby01));

  // relativistic BW + Cheby
  RooAddPdf kpmass_qq_pdf("kpmass_qq_pdf", "kpmass qq pdf", RooArgList(kpmass_qq_bw_pdf, kpmass_qq_cheby_pdf), bw_qq_fraction2);


  //----------------------------------- k0 hcos ---------------------------------------------
    
  //----------------------
  // 0.78 < m(K*0) < 0.83
  //----------------------

  // Chebychev

  RooRealVar k0hcos_qq_c0( "k0hcos_qq_c0","a_{0}(Cheby)", 0.656);
  RooRealVar k0hcos_qq_c1( "k0hcos_qq_c1","a_{1}(Cheby)",-0.1565);
  RooRealVar k0hcos_qq_c2( "k0hcos_qq_c2","a_{2}(Cheby)", 0.171);
  RooRealVar k0hcos_qq_c3( "k0hcos_qq_c3","a_{3}(Cheby)",-0.0053);
  RooRealVar k0hcos_qq_c4( "k0hcos_qq_c4","a_{4}(Cheby)",-0.0050);
  RooRealVar k0hcos_qq_c5( "k0hcos_qq_c5","a_{5}(Cheby)", 0.145);

  // error function
  RooRealVar k0hcos_qq_interxc( "k0hcos_qq_interxc","x intercept", 0.9797);
  RooRealVar k0hcos_qq_sigmac( "k0hcos_qq_sigmac","sigma", 0.0366);

  RooFormulaVar k0hcos_qq_effc     ("k0hcos_qq_effc",      "0.5*(TMath::Erf((-k0_hcos+k0hcos_qq_interxc)/k0hcos_qq_sigmac)+1.0)",RooArgSet(k0_hcos, k0hcos_qq_interxc, k0hcos_qq_sigmac));
  RooChebychev  k0hcos_qq_chev_pdfc("k0hcos_qq_chev_pdfc", "pdf", k0_hcos, RooArgList(k0hcos_qq_c0, k0hcos_qq_c1, k0hcos_qq_c2, k0hcos_qq_c3, k0hcos_qq_c4, k0hcos_qq_c5));
  RooEffProd    k0hcos_qq_pdfc     ("k0hcos_qq_pdfc",      "model with efficiency", k0hcos_qq_chev_pdfc, k0hcos_qq_effc) ;
  

  //----------------------
  // 0.83 < m(K*0) < 1.00
  //----------------------

  RooRealVar k0hcos_qq_d0( "k0hcos_qq_d0","a_{0}(Cheby)", 0.1316);
  RooRealVar k0hcos_qq_d1( "k0hcos_qq_d1","a_{1}(Cheby)", 0.0611);
  RooRealVar k0hcos_qq_d2( "k0hcos_qq_d2","a_{2}(Cheby)",-0.07371);
  RooRealVar k0hcos_qq_d3( "k0hcos_qq_d3","a_{3}(Cheby)",-0.05088);

  RooRealVar k0hcos_qq_interxd( "k0hcos_qq_interxd","x intercept", 0.95893);
  RooRealVar k0hcos_qq_sigmad( "k0hcos_qq_sigmad","sigma", 0.01585);

  RooFormulaVar k0hcos_qq_effd     ("k0hcos_qq_effd",      "0.5*(TMath::Erf((-k0_hcos+k0hcos_qq_interxd)/k0hcos_qq_sigmad)+1.0)",RooArgSet(k0_hcos, k0hcos_qq_interxd, k0hcos_qq_sigmad));
  RooChebychev  k0hcos_qq_chev_pdfd("k0hcos_qq_chev_pdfd", "pdf", k0_hcos, RooArgList(k0hcos_qq_d0, k0hcos_qq_d1, k0hcos_qq_d2, k0hcos_qq_d3));
  RooEffProd    k0hcos_qq_pdfd     ("k0hcos_qq_pdfd",      "model with efficiency", k0hcos_qq_chev_pdfd, k0hcos_qq_effd) ;

  //
  // Simultaneous Fit
  //
  RooSimultaneous k0hcos_qq_pdf("k0hcos_qq_pdf","simultaneous pdf",sample) ;
  k0hcos_qq_pdf.addPdf(k0hcos_qq_pdfc, "r01");
  k0hcos_qq_pdf.addPdf(k0hcos_qq_pdfd, "r02");


  //----------------------------------- kp hcos ---------------------------------------------

  // Cheby
  RooRealVar kphcos_qq_c0( "kphcos_qq_c0","a_{0}(Cheby)",-0.18729);
  RooRealVar kphcos_qq_c1( "kphcos_qq_c1","a_{1}(Cheby)", 0.0649);
  RooRealVar kphcos_qq_c2( "kphcos_qq_c2","a_{2}(Cheby)", 0.0617);
  RooRealVar kphcos_qq_c3( "kphcos_qq_c3","a_{3}(Cheby)",-0.03304);
  RooRealVar kphcos_qq_c4( "kphcos_qq_c4","a_{4}(Cheby)", 0.0262);

  // error function
  RooRealVar kphcos_qq_interx( "kphcos_qq_interx","x intercept", 0.96145);
  RooRealVar kphcos_qq_sigma ( "kphcos_qq_sigma","sigma", 0.0265);

  RooFormulaVar kphcos_qq_eff("kphcos_qq_eff","0.5*(TMath::Erf((kp_hcos+kphcos_qq_interx)/kphcos_qq_sigma)+1.0)", RooArgSet(kp_hcos,kphcos_qq_interx,kphcos_qq_sigma) );
  RooChebychev kphcos_qq_chev_pdf("kphcos_qq_chev_pdf", "pdf", kp_hcos, RooArgList(kphcos_qq_c0,kphcos_qq_c1,kphcos_qq_c2,kphcos_qq_c3,kphcos_qq_c4));
  RooEffProd kphcos_qq_pdf("kphcos_qq_pdf","model with efficiency",kphcos_qq_chev_pdf,kphcos_qq_eff) ;


  //----------------------------------- NN out ----------------------------------------------

  // Gaussian1
  RooRealVar    nn_tran2_qq_gau1_mean("mean_{G1}",  "mean(NB_{trans},Gau1)",  -0.4456, -2.5, 1.5);
  RooRealVar    nn_tran2_qq_gau1_wid ("width_{G1}", "#sigma(NB_{trans},Gau1)", 1.5499,  0.0, 3.0);

  // Gaussian2 
  //  RooRealVar    nn_tran2_qq_gau2_mean("mean_{G2}",  "mean(NB_{trans},Gau1)",  -2.3170, -3.5,-1.5);
  //  RooRealVar    nn_tran2_qq_gau2_wid ("width_{G2}", "#sigma(NB_{trans},Gau1)", 2.0590,  1.0, 3.0);
  RooFormulaVar nn_tran2_qq_gau2_mean("nn_tran2_qq_gau2_mean", "mean_{G1}-1.8714",  nn_tran2_qq_gau1_mean);
  RooFormulaVar nn_tran2_qq_gau2_wid ("nn_tran2_qq_gau2_wid",  "width_{G1}*1.3285", nn_tran2_qq_gau1_wid);

  // fraction between Gaussian1 and Gaussian2
  RooRealVar nn_tran2_qq_frac("nn_tran2_qq_frac", "nn_tran2 qq fraction 2 Bifurs", 0.8630, 0.7, 1.0);

  RooGaussian nn_tran2_qq_gau1_pdf("nn_tran2_qq_gau1_pdf", "nn_tran2 qq Gaussian", nn_tran2, nn_tran2_qq_gau1_mean, nn_tran2_qq_gau1_wid);
  RooGaussian nn_tran2_qq_gau2_pdf("nn_tran2_qq_gau2_pdf", "nn_tran2 qq Gaussian", nn_tran2, nn_tran2_qq_gau2_mean, nn_tran2_qq_gau2_wid);

  // Gaussian1 and Gaussian2
  RooAddPdf nn_tran2_qq_pdf ("nn_tran2_qq_pdf", "nn_tran2 qq pdf", RooArgList(nn_tran2_qq_gau1_pdf, nn_tran2_qq_gau2_pdf), nn_tran2_qq_frac);

  //
  // bb PDF Setting
  //

  //-------------------------------------- de -----------------------------------------------

  // Cheby
  RooRealVar   de_bb_cheby0( "de_bb_cheby0","a_{0}(#DeltaE,Cheby)", -0.5618);

  RooChebychev de_bb_pdf("de_bb_pdf", "de bb pdf", de, RooArgList(de_bb_cheby0));


  //----------------------------------- mbc ----------------------------------------------

  // Argus - mbc(bb)
  RooRealVar mbc_bb_argus_par( "mbc_bb_argus_par", "mbc bb argus parameter", -43.27);

  RooArgusBG mbc_bb_pdf( "mbc_bb_pdf", "mbc bb pdf", mbc, mbc_qq_argus_cutoff, mbc_bb_argus_par);


  //----------------------------------- k0 mass ---------------------------------------------

  // Cheby
  RooRealVar   k0mass_bb_cheby( "k0mass_bb_cheby","a_{0}(m(K*0),Cheby)", 0.371);

  RooChebychev k0mass_bb_pdf("k0mass_bb_pdf", "k0mass bb pdf", k0_mass, RooArgList(k0mass_bb_cheby));

  //----------------------------------- kp mass ---------------------------------------------

  // Cheby 
  RooRealVar   kpmass_bb_cheby0( "kpmass_bb_cheby0","a_{0}(m(K*+),Cheby)", 0.125);
  RooRealVar   kpmass_bb_cheby1( "kpmass_bb_cheby1","a_{0}(m(K*+),Cheby)",-0.0878);

  RooChebychev kpmass_bb_pdf("kpmass_bb_pdf", "kpmass bb pdf", kp_mass, RooArgList(kpmass_bb_cheby0,kpmass_bb_cheby1));


  //----------------------------------- k0 hcos ---------------------------------------------
  
  //----------------------
  // 0.78 < m(K*0) < 0.83
  //----------------------

  // Cheby 
  RooRealVar k0hcos_bb_c0( "k0hcos_bb_c0","a_{0}(Cheby)", 0.954);
  RooRealVar k0hcos_bb_c1( "k0hcos_bb_c1","a_{1}(Cheby)", 0.182);
  RooRealVar k0hcos_bb_c2( "k0hcos_bb_c2","a_{2}(Cheby)",-0.0085);
  RooRealVar k0hcos_bb_c3( "k0hcos_bb_c3","a_{3}(Cheby)",-0.2281);

  RooChebychev  k0hcos_bb_pdfc("k0hcos_bb_pdfc", "pdf", k0_hcos, RooArgList(k0hcos_bb_c0,k0hcos_bb_c1,k0hcos_bb_c2,k0hcos_bb_c3));

  //----------------------
  // 0.83 < m(K*0) < 1.00
  //----------------------

  // Cheby 
  RooRealVar k0hcos_bb_d0( "k0hcos_bb_d0","a_{0}(Cheby)", 0.252);
  RooRealVar k0hcos_bb_d1( "k0hcos_bb_d1","a_{1}(Cheby)", 0.685);
  RooRealVar k0hcos_bb_d2( "k0hcos_bb_d2","a_{2}(Cheby)", 0.007);
  RooRealVar k0hcos_bb_d3( "k0hcos_bb_d3","a_{3}(Cheby)", 0.070);

  // error function
  RooRealVar k0hcos_bb_interx( "k0hcos_bb_interx","x intercept", 0.930);
  RooRealVar k0hcos_bb_sigma ( "k0hcos_bb_sigma", "sigma",       0.037);

  RooFormulaVar k0hcos_bb_eff     ("k0hcos_bb_eff",      "0.5*(TMath::Erf((-k0_hcos+k0hcos_bb_interx)/k0hcos_bb_sigma)+1.0)",RooArgSet(k0_hcos,k0hcos_bb_interx,k0hcos_bb_sigma));
  RooChebychev  k0hcos_bb_chev_pdf("k0hcos_bb_chev_pdf", "pdf", k0_hcos, RooArgList(k0hcos_bb_d0,k0hcos_bb_d1,k0hcos_bb_d2,k0hcos_bb_d3));
  RooEffProd    k0hcos_bb_pdfd    ("k0hcos_bb_pdfd",     "model with efficiency",k0hcos_bb_chev_pdf,k0hcos_bb_eff) ;
  

  //
  // Simultaneous Fit
  //
  RooSimultaneous k0hcos_bb_pdf("k0hcos_bb_pdf","simultaneous pdf",sample) ;
  k0hcos_bb_pdf.addPdf(k0hcos_bb_pdfc, "r01");
  k0hcos_bb_pdf.addPdf(k0hcos_bb_pdfd, "r02");


  //----------------------------------- kp hcos ---------------------------------------------

  // Cheby 
  RooRealVar kphcos_bb_c0( "kphcos_bb_c0","a_{0}(Cheby)",-0.0586);
  RooRealVar kphcos_bb_c1( "kphcos_bb_c1","a_{1}(Cheby)",-0.0091);
  RooRealVar kphcos_bb_c2( "kphcos_bb_c2","a_{2}(Cheby)", 0.2980);
  RooRealVar kphcos_bb_c3( "kphcos_bb_c3","a_{3}(Cheby)",-0.1854);
  RooRealVar kphcos_bb_c4( "kphcos_bb_c4","a_{4}(Cheby)", 0.1560);

  RooChebychev kphcos_bb_pdf("kphcos_bb_pdf", "pdf", kp_hcos, RooArgList(kphcos_bb_c0,kphcos_bb_c1,kphcos_bb_c2,kphcos_bb_c3,kphcos_bb_c4));

  //----------------------------------- NN out ----------------------------------------------

  // BifurGauss 
  RooRealVar nn_tran2_bb_bif_mean("nn_tran2_bb_bif_mean", "nn_tran2 bb mean",    1.79);
  RooRealVar nn_tran2_bb_bif_widL("nn_tran2_bb_bif_widL", "nn_tran2 bb width L", 2.435);
  RooRealVar nn_tran2_bb_bif_widR("nn_tran2_bb_bif_widR", "nn_tran2 bb width R", 1.365);

  RooBifurGauss nn_tran2_bb_pdf("nn_tran2_bb_pdf",  "nn_tran2 bb BifurGauss", nn_tran2, nn_tran2_bb_bif_mean, nn_tran2_bb_bif_widL, nn_tran2_bb_bif_widR);


  //
  // rare PDF Setting
  //

  //-------------------------------------- de -----------------------------------------------


  RooRealVar de_rare_cheby0("de_rare_cheby0","a0", -0.5027);  
  RooRealVar de_rare_cheby1("de_rare_cheby1","a1",  0.1520);
  RooRealVar de_rare_fR("de_rare_fR",  "fr", 0.101);
  RooRealVar de_rare_gauR_mean("de_rare_gauR_mean", "mean",   0.0653);
  RooRealVar de_rare_gauR_wid ("de_rare_gauR_wid",  "width",  0.0254);

  RooGaussian  de_rare_pdfR("de_rare_pdfR", "Gaussian", de, de_rare_gauR_mean, de_rare_gauR_wid);
  RooChebychev de_rare_pdf0("de_rare_pdf0", "BKG",      de, RooArgList(de_rare_cheby0,de_rare_cheby1));

  RooAddPdf de_rare_pdf("de_rare_pdf", "pdf", RooArgList(de_rare_pdfR, de_rare_pdf0), de_rare_fR);


  //----------------------------------- mbc ----------------------------------------------

  // Gaussian
  RooRealVar mbc_rare_mean  ("mbc_rare_mean",   "mean",        5.27917);
  RooRealVar mbc_rare_wid   ("mbc_rare_wid",    "width",       0.00303);

  // Argus
  RooRealVar mbc_rare_par   ("mbc_rare_par",    "parameter", -64.42);
  RooRealVar mbc_rare_fg    ("mbc_rare_fg",     "fraction",    0.209);

  // Gaussian + Argus  
  RooGaussian mbc_rare_pdf0("mbc_rare_pdf0",  "true signal Gaussian", mbc, mbc_rare_mean, mbc_rare_wid);
  RooArgusBG  mbc_rare_bkg( "mbc_rare_bkg",   "bkg",                  mbc, mbc_qq_argus_cutoff, mbc_rare_par);

  RooAddPdf mbc_rare_pdf("mbc_rare_pdf", "pdf", RooArgList(mbc_rare_pdf0, mbc_rare_bkg), mbc_rare_fg);


  //----------------------------------- k0 mass ---------------------------------------------

  // BW
  RooRealVar k0mass_rare_bw_mean("k0mass_rare_bw_mean", "k0mass rare mean",   0.8900);
  RooRealVar k0mass_rare_bw_wid ("k0mass_rare_bw_wid",  "k0mass rare width",  0.0540);

  // Cheby
  RooRealVar k0mass_rare_cheby01("k0mass_rare_cheby01", "a_{0}(m(K*0),Cheby)", 0.298);

  // fraction between BW+Cheby
  RooRealVar bw_rare_fraction1  ("bw_rare_fraction1",   "bw rare fraction",    0.231);

  RooBreitWigner k0mass_rare_bw_pdf   ("k0mass_rare_bw_pdf",    "k0mass rare BW",         k0_mass, k0mass_rare_bw_mean, k0mass_rare_bw_wid);  
  RooChebychev   k0mass_rare_cheby_pdf("k0mass_rare_cheby_pdf", "k0_mass rare cheby pdf", k0_mass, RooArgList(k0mass_rare_cheby01));

  // BW+Cheby
  RooAddPdf k0mass_rare_pdf("k0mass_rare_pdf", "k0mass rare pdf", RooArgList(k0mass_rare_bw_pdf, k0mass_rare_cheby_pdf), bw_rare_fraction1);


  //----------------------------------- kp mass ---------------------------------------------

  // BW
  RooRealVar kpmass_rare_bw_mean("kpmass_rare_bw_mean", "kpmass rare mean",   0.8939);
  RooRealVar kpmass_rare_bw_wid ("kpmass_rare_bw_wid",  "kpmass rare width",  0.0533);

  // Cheby
  RooRealVar kpmass_rare_cheby01("kpmass_rare_cheby01", "a_{0}(m(K*+),Cheby)", 0.088);

  // fraction between BW+Cheby
  RooRealVar bw_rare_fraction2  ("bw_rare_fraction2",   "bw rare fraction",    0.295);

  RooBreitWigner kpmass_rare_bw_pdf   ("kpmass_rare_bw_pdf",    "kpmass rare BW",         kp_mass, kpmass_rare_bw_mean, kpmass_rare_bw_wid);  
  RooChebychev   kpmass_rare_cheby_pdf("kpmass_rare_cheby_pdf", "kp_mass rare cheby pdf", kp_mass, RooArgList(kpmass_rare_cheby01));

  // BW+Cheby
  RooAddPdf kpmass_rare_pdf("kpmass_rare_pdf", "kpmass rare pdf", RooArgList(kpmass_rare_bw_pdf, kpmass_rare_cheby_pdf), bw_rare_fraction2);


  //----------------------------------- k0 hcos ---------------------------------------------

                               /*   HistPdf modeling   */

  RooSimultaneous k0hcos_rare_pdf("k0hcos_rare_pdf","simultaneous pdf",sample) ;
  k0hcos_rare_pdf.addPdf(*model59, "r01");
  k0hcos_rare_pdf.addPdf(*model60, "r02");

  //----------------------------------- kp hcos ---------------------------------------------

                               /*   HistPdf modeling   */


  //----------------------------------- NN out ----------------------------------------------

  // BifurGauss 
  RooRealVar nn_tran2_rare_bif_mean("nn_tran2_rare_bif_mean", "nn_tran2 rare mean",    2.366);
  RooRealVar nn_tran2_rare_bif_widL("nn_tran2_rare_bif_widL", "nn_tran2 rare width L", 2.558);
  RooRealVar nn_tran2_rare_bif_widR("nn_tran2_rare_bif_widR", "nn_tran2 rare width R", 1.400);

  RooBifurGauss nn_tran2_rare_pdf("nn_tran2_rare_pdf", "nn_tran2 rare BifurGauss", nn_tran2, nn_tran2_rare_bif_mean, nn_tran2_rare_bif_widL, nn_tran2_rare_bif_widR);





  //
  //  added pdf 
  //
  RooRealVar nsig ("nsig", "N_{sig}", vnsig,  -150.0,   850.0);
  RooRealVar fl   ("fl",   "f_{L}",   ifl,      -4.5,     6.5);
  RooRealVar nqq  ("nqq",  "N_{qq}",  vnqq,  17000.0, 22000.0);
  RooRealVar nbb  ("nbb",  "N_{bb}",  vnbb,   -200.0,   700.0);
  nbb.setConstant(0);
  RooRealVar nrare("nrare","N_{rare}",vnrare);

  RooFormulaVar nsig0("nsig0", "nsig*(1.0-fl)", RooArgList(nsig,fl));
  RooFormulaVar nsig1("nsig1", "nsig*fl",       RooArgList(nsig,fl));

  //
  //  product pdf
  //
  RooProdPdf final_true_sig_pdf0("final_true_sig_pdf0", "final tr true sig pdf", RooArgList(de_true_sig_pdf0, mbc_true_sig_pdf0, k0mass_true_sig_pdf0, kpmass_true_sig_pdf0, k0hcos_true_sig_pdf0, *model62, nn_tran2_true_sig_pdf0));
  RooProdPdf final_true_sig_pdf1("final_true_sig_pdf1", "final lo true sig pdf", RooArgList(de_true_sig_pdf1, mbc_true_sig_pdf1, k0mass_true_sig_pdf1, kpmass_true_sig_pdf1, k0hcos_true_sig_pdf1, *model61, nn_tran2_true_sig_pdf1));

  RooProdPdf final_scf_sig_pdf0("final_scf_sig_pdf0", "final tr scf sig pdf", RooArgList( *model12, *model22,  *model32, *model42, k0hcos_scf_sig_pdf0, *model64, nn_tran2_scf_sig_pdf0));
  RooProdPdf final_scf_sig_pdf1("final_scf_sig_pdf1", "final lo scf sig pdf", RooArgList( *model11, *model21,  *model31, *model41, k0hcos_scf_sig_pdf1, *model63, nn_tran2_scf_sig_pdf1));

  RooAddPdf  final_sig_pdf0 ("final_sig_pdf0",  "final tr signal pdf",  RooArgList(final_true_sig_pdf0, final_scf_sig_pdf0), true_fraction0 );
  RooAddPdf  final_sig_pdf1 ("final_sig_pdf1",  "final lo signal pdf",  RooArgList(final_true_sig_pdf1, final_scf_sig_pdf1), true_fraction1 );

  RooAddPdf  final_sig_pdf ("final_sig_pdf",  "final signalpdf",  RooArgList(final_sig_pdf0, final_sig_pdf1), RooArgList(nsig0,nsig1) );


  RooProdPdf final_qq_pdf  ("final_qq_pdf",   "final qq pdf",     RooArgList(de_qq_pdf,   mbc_qq_pdf,   k0mass_qq_pdf,   kpmass_qq_pdf,   k0hcos_qq_pdf,   kphcos_qq_pdf,   nn_tran2_qq_pdf)   );
  RooProdPdf final_bb_pdf  ("final_bb_pdf",   "final bb pdf",     RooArgList(de_bb_pdf,   mbc_bb_pdf,   k0mass_bb_pdf,   kpmass_bb_pdf,   k0hcos_bb_pdf,   kphcos_bb_pdf,   nn_tran2_bb_pdf)   );
  RooProdPdf final_rare_pdf("final_rare_pdf", "final rare pdf",   RooArgList(de_rare_pdf, mbc_rare_pdf, k0mass_rare_pdf, kpmass_rare_pdf, k0hcos_rare_pdf, *model67,        nn_tran2_rare_pdf) );
  //
  //
  //
  RooProdPdf final_true_sig_pdf0a("final_true_sig_pdf0a", "final true tr sig pdf", RooArgList(de_true_sig_pdf0, mbc_true_sig_pdf0, k0mass_true_sig_pdf0, kpmass_true_sig_pdf0, *model53, *model62, nn_tran2_true_sig_pdf0));
  RooProdPdf final_true_sig_pdf0b("final_true_sig_pdf0b", "final true tr sig pdf", RooArgList(de_true_sig_pdf0, mbc_true_sig_pdf0, k0mass_true_sig_pdf0, kpmass_true_sig_pdf0, *model54, *model62, nn_tran2_true_sig_pdf0));
  RooProdPdf final_scf_sig_pdf0a("final_scf_sig_pdf0a", "final scf tr sig pdf", RooArgList(*model12, *model22, *model32, *model42, *model57, *model64, nn_tran2_scf_sig_pdf0));
  RooProdPdf final_scf_sig_pdf0b("final_scf_sig_pdf0b", "final scf tr sig pdf", RooArgList(*model12, *model22, *model32, *model42, *model58, *model64, nn_tran2_scf_sig_pdf0));

  RooProdPdf final_true_sig_pdf1a("final_true_sig_pdf1a", "final true lo sig pdf", RooArgList(de_true_sig_pdf1, mbc_true_sig_pdf1, k0mass_true_sig_pdf1, kpmass_true_sig_pdf1, *model51, *model61, nn_tran2_true_sig_pdf1));
  RooProdPdf final_true_sig_pdf1b("final_true_sig_pdf1b", "final true lo sig pdf", RooArgList(de_true_sig_pdf1, mbc_true_sig_pdf1, k0mass_true_sig_pdf1, kpmass_true_sig_pdf1, *model52, *model61, nn_tran2_true_sig_pdf1));
  RooProdPdf final_scf_sig_pdf1a("final_scf_sig_pdf1a", "final scf lo sig pdf", RooArgList(*model11, *model21, *model31, *model41, *model55, *model63, nn_tran2_scf_sig_pdf1));
  RooProdPdf final_scf_sig_pdf1b("final_scf_sig_pdf1b", "final scf lo sig pdf", RooArgList(*model11, *model21, *model31, *model41, *model56, *model63, nn_tran2_scf_sig_pdf1));

  RooAddPdf  final_sig_pdf0a ("final_sig_pdf0a",  "final tr signal pdf",  RooArgList(final_true_sig_pdf0a, final_scf_sig_pdf0a), true_fraction0A );
  RooAddPdf  final_sig_pdf0b ("final_sig_pdf0b",  "final tr signal pdf",  RooArgList(final_true_sig_pdf0b, final_scf_sig_pdf0b), true_fraction0B );
  RooAddPdf  final_sig_pdf1a ("final_sig_pdf1a",  "final lo signal pdf",  RooArgList(final_true_sig_pdf1a, final_scf_sig_pdf1a), true_fraction1A );
  RooAddPdf  final_sig_pdf1b ("final_sig_pdf1b",  "final lo signal pdf",  RooArgList(final_true_sig_pdf1b, final_scf_sig_pdf1b), true_fraction1B );

  RooAddPdf  final_sig_pdfa ("final_sig_pdfa",  "final signalpdf",  RooArgList(final_sig_pdf0a, final_sig_pdf1a), RooArgList(nsig0,nsig1) );
  RooAddPdf  final_sig_pdfb ("final_sig_pdfb",  "final signalpdf",  RooArgList(final_sig_pdf0b, final_sig_pdf1b), RooArgList(nsig0,nsig1) );

  RooProdPdf final_qq_pdfa ("final_qq_pdfa",   "final qq pdf",     RooArgList(de_qq_pdf,   mbc_qq_pdf,   k0mass_qq_pdf,   kpmass_qq_pdf,   k0hcos_qq_pdfc,  kphcos_qq_pdf,   nn_tran2_qq_pdf)   );
  RooProdPdf final_qq_pdfb ("final_qq_pdfb",   "final qq pdf",     RooArgList(de_qq_pdf,   mbc_qq_pdf,   k0mass_qq_pdf,   kpmass_qq_pdf,   k0hcos_qq_pdfd,  kphcos_qq_pdf,   nn_tran2_qq_pdf)   );

  RooProdPdf final_bb_pdfa ("final_bb_pdfa",   "final bb pdf",     RooArgList(de_bb_pdf,   mbc_bb_pdf,   k0mass_bb_pdf,   kpmass_bb_pdf,   k0hcos_bb_pdfc,  kphcos_bb_pdf,   nn_tran2_bb_pdf)   );
  RooProdPdf final_bb_pdfb ("final_bb_pdfb",   "final bb pdf",     RooArgList(de_bb_pdf,   mbc_bb_pdf,   k0mass_bb_pdf,   kpmass_bb_pdf,   k0hcos_bb_pdfd,  kphcos_bb_pdf,   nn_tran2_bb_pdf)   );

  RooProdPdf final_rare_pdfa("final_rare_pdfa", "final rare pdf",  RooArgList(de_rare_pdf, mbc_rare_pdf, k0mass_rare_pdf, kpmass_rare_pdf, *model59,        *model67,        nn_tran2_rare_pdf) );
  RooProdPdf final_rare_pdfb("final_rare_pdfb", "final rare pdf",  RooArgList(de_rare_pdf, mbc_rare_pdf, k0mass_rare_pdf, kpmass_rare_pdf, *model60,        *model67,        nn_tran2_rare_pdf) );

  //
  //  selected pdf list
  //
  RooArgList* pdflist;
  RooArgList* parlist;
  cout<<" "<<endl;
  cout<<" "<<endl;
  if( set_sig>0 && set_qq>0 && set_bb==0 &&set_rare==0 ) {
    pdflist = new RooArgList(final_sig_pdf0, final_sig_pdf1, final_qq_pdf); 
    parlist = new RooArgList(         nsig0,          nsig1,          nqq); 
    cout<<"p.d.f. list and yield list :  signal + qq"<<endl;
  }
  if( set_sig>0 && set_bb>0 && set_qq==0 &&set_rare==0   ) { 
    pdflist = new RooArgList(final_sig_pdf0, final_sig_pdf1, final_bb_pdf); 
    parlist = new RooArgList(         nsig0,          nsig1,          nbb); 
    cout<<"p.d.f. list and yield list :  signal + bb"<<endl;
  }
  if( set_sig>0 && set_rare>0 && set_qq==0 &&set_bb==0 ) { 
    pdflist = new RooArgList(final_sig_pdf0, final_sig_pdf1, final_rare_pdf); 
    parlist = new RooArgList(         nsig0,          nsig1,        nrare); 
    cout<<"p.d.f. list and yield list :  signal + rare"<<endl;
  }
  if( set_qq>0  && set_bb>0  && set_sig==0 &&set_rare==0  ) { 
    pdflist = new RooArgList(final_qq_pdf, final_bb_pdf); 
    parlist = new RooArgList(         nqq,          nbb); 
    cout<<"p.d.f. list and yield list :  bb + qq"<<endl;
  }
  if( set_qq>0  && set_rare>0 && set_sig==0 &&set_bb==0 ) { 
    pdflist = new RooArgList(final_qq_pdf, final_rare_pdf); 
    parlist = new RooArgList(         nqq,          nrare); 
    cout<<"p.d.f. list and yield list :  qq + rare"<<endl;
  }
  if( set_bb>0  && set_rare>0 && set_sig==0 &&set_qq==0 ) { 
    pdflist = new RooArgList(final_bb_pdf, final_rare_pdf); 
    parlist = new RooArgList(         nbb,          nrare); 
    cout<<"p.d.f. list and yield list :  bb + rare"<<endl;
  }
  if( set_sig>0 && set_qq>0 && set_bb>0 && set_rare==0   ) { 
    pdflist = new RooArgList(final_sig_pdf0, final_sig_pdf1, final_qq_pdf, final_bb_pdf); 
    parlist = new RooArgList(         nsig0,          nsig1,          nqq,          nbb); 
    cout<<"p.d.f. list and yield list :  signal + qq + bb"<<endl;
  }
  if( set_sig>0 && set_qq>0 && set_rare>0 && set_bb==0 ) { 
    pdflist = new RooArgList(final_sig_pdf0, final_sig_pdf1, final_qq_pdf, final_rare_pdf); 
    parlist = new RooArgList(         nsig0,          nsig1,          nqq,          nrare); 
    cout<<"p.d.f. list and yield list :  signal + qq + rare"<<endl;
  }
  if( set_sig>0 && set_bb>0 && set_rare>0 && set_qq==0 ) { 
    pdflist = new RooArgList(final_sig_pdf0, final_sig_pdf1, final_bb_pdf, final_rare_pdf); 
    parlist = new RooArgList(         nsig0,          nsig1,          nbb,          nrare); 
    cout<<"p.d.f. list and yield list :  signal + bb + rare"<<endl;
  }
  if( set_qq>0  && set_bb>0 && set_rare>0 && set_sig==0 ) { 
    pdflist = new RooArgList(final_qq_pdf,   final_bb_pdf,   final_rare_pdf); 
    parlist = new RooArgList(         nqq,            nbb,            nrare); 
    cout<<"p.d.f. list and yield list :  qq + bb + rare"<<endl;
  }
  if( set_sig>0 && set_qq>0 && set_bb>0 && set_rare>0 ) { 
    pdflist = new RooArgList(final_sig_pdf0, final_sig_pdf1, final_qq_pdf, final_bb_pdf, final_rare_pdf); 
    parlist = new RooArgList(         nsig0,          nsig1,          nqq,          nbb,          nrare); 
    cout<<"p.d.f. list and yield list :  signal + qq + bb + rare"<<endl;
  }
  //
  cout<<" "<<endl;


  //
  // final pdf
  //
  RooAddPdf final_pdf("final_pdf", "final pdf", *pdflist, *parlist );

  //
  // close opened files
  //
  f11->Close("R");
  f12->Close("R");
  f21->Close("R");
  f22->Close("R");
  f31->Close("R");
  f32->Close("R");
  f41->Close("R");
  f42->Close("R");
  f51->Close("R");
  f52->Close("R");
  f53->Close("R");
  f54->Close("R");
  f55->Close("R");
  f56->Close("R");
  f57->Close("R");
  f58->Close("R");
  f59->Close("R");
  f60->Close("R");
  f61->Close("R");
  f62->Close("R");
  f63->Close("R");
  f64->Close("R");
  f67->Close("R");

  //
  // delete memories of pointer
  //
  delete f11;
  delete f12;
  delete f21;
  delete f22;
  delete f31;
  delete f32;
  delete f41;
  delete f42;
  delete f51;
  delete f52;
  delete f53;
  delete f54;
  delete f55;
  delete f56;
  delete f57;
  delete f58;
  delete f59;
  delete f60;
  delete f61;
  delete f62;
  delete f63;
  delete f64;
  delete f67;


  //
  // DataSet 
  //
  TChain* xxx = new TChain("h1"); 
  if(set_sig==1) {  xxx->Add(in_root1.c_str());  cout<<" "<<endl; cout<<"GSim signal sample is involved!"<<endl;   cout<<" "<<endl; }
  if(set_qq==1)  {  xxx->Add(in_root2.c_str());  cout<<" "<<endl; cout<<"GSim qq     sample is involved!"<<endl;   cout<<" "<<endl; }
  if(set_bb==1)  {  xxx->Add(in_root3.c_str());  cout<<" "<<endl; cout<<"GSim bb     sample is involved!"<<endl;   cout<<" "<<endl; }
  if(set_rare==1){  xxx->Add(in_root4.c_str());  cout<<" "<<endl; cout<<"GSim rare   sample is involved!"<<endl;   cout<<" "<<endl; }
  xxx->Add("zeroevent.root");

  TTree* yyy01;
  TTree* yyy02;
  TChain* yyy03 = new TChain("h1"); 
  TChain* yyy04 = new TChain("h1"); 
  RooDataSet* data01;
  RooDataSet* data02;
  if(set_sig==1||set_qq==1||set_bb==1||set_rare==1){
    yyy01 = xxx->CopyTree("k0_mass>0.78&& k0_mass<0.83");
    yyy02 = xxx->CopyTree("k0_mass>0.83&& k0_mass<1.00");
    data01 = new RooDataSet("data01", "Data", yyy01, *varset);
    data02 = new RooDataSet("data02", "Data", yyy02, *varset);
  }else{
    yyy03->Add("zeroevent.root");
    yyy04->Add("zeroevent.root");
    data01 = new RooDataSet("data01", "Data", yyy03, *varset);
    data02 = new RooDataSet("data02", "Data", yyy04, *varset);
  }

  //
  // ToyMC dataset & generation
  //
  TRandom3* rand = new TRandom3(0);

  int tmp_nsig  = rand->Poisson(vnsig);
  int tmp_nqq   = rand->Poisson(vnqq);
  int tmp_nbb   = rand->Poisson(vnbb);
  int tmp_nrare = rand->Poisson(vnrare);

  RooDataSet* data_siga  = final_sig_pdfa.generate(*varset,tmp_nsig);
  RooDataSet* data_sigb  = final_sig_pdfb.generate(*varset,tmp_nsig);
  RooDataSet* data_qqa   = final_qq_pdfa.generate(*varset,tmp_nqq);
  RooDataSet* data_qqb   = final_qq_pdfb.generate(*varset,tmp_nqq);
  RooDataSet* data_bba   = final_bb_pdfa.generate(*varset,tmp_nbb);
  RooDataSet* data_bbb   = final_bb_pdfb.generate(*varset,tmp_nbb);
  RooDataSet* data_rarea = final_rare_pdfa.generate(*varset,tmp_nrare);
  RooDataSet* data_rareb = final_rare_pdfb.generate(*varset,tmp_nrare);

  cout<<" "<<endl;
  cout<<" "<<endl;

  delete rand;

  RooDataSet* rdc_data_sig01  = (RooDataSet*) data_siga->reduce( "k0_mass>0.78&& k0_mass<0.83");
  RooDataSet* rdc_data_qq01   = (RooDataSet*) data_qqa->reduce(  "k0_mass>0.78&& k0_mass<0.83");
  RooDataSet* rdc_data_bb01   = (RooDataSet*) data_bba->reduce(  "k0_mass>0.78&& k0_mass<0.83");
  RooDataSet* rdc_data_rare01 = (RooDataSet*) data_rarea->reduce("k0_mass>0.78&& k0_mass<0.83");
  RooDataSet* rdc_data_sig02  = (RooDataSet*) data_sigb->reduce( "k0_mass>0.83&& k0_mass<1.00");
  RooDataSet* rdc_data_qq02   = (RooDataSet*) data_qqb->reduce(  "k0_mass>0.83&& k0_mass<1.00");
  RooDataSet* rdc_data_bb02   = (RooDataSet*) data_bbb->reduce(  "k0_mass>0.83&& k0_mass<1.00");
  RooDataSet* rdc_data_rare02 = (RooDataSet*) data_rareb->reduce("k0_mass>0.83&& k0_mass<1.00");


  if(set_sig==2) {  
    data01->append(*rdc_data_sig01);
    data02->append(*rdc_data_sig02);
    cout<<"Toy signal sample is generated!"<<endl;  
  }
  if(set_qq==2)  {  
    data01->append(*rdc_data_qq01);
    data02->append(*rdc_data_qq02);
    cout<<"Toy qq     sample is generated!"<<endl; 
  }
  if(set_bb==2)  {  
    data01->append(*rdc_data_bb01);
    data02->append(*rdc_data_bb02);
    cout<<"Toy bb     sample is generated!"<<endl;  
  }
  if(set_rare==2){  
    data01->append(*rdc_data_rare01);
    data02->append(*rdc_data_rare02);
    cout<<"Toy rare   sample is generated!"<<endl;  
  }


  RooDataSet combData("combData", "combined data", *varset, Index(sample), Import("r01", *data01), Import("r02", *data02) ) ;


  //
  // Fit
  //
  RooFitResult *final_result = new RooFitResult("final_result", "final result");
  final_result = final_pdf.fitTo( combData, Hesse(1), Minos(0), Timer(0), Extended(true), Save(0));


  //
  // Save results
  //
  float v_nsig                  = nsig.getVal();
  float e_nsig                  = nsig.getError(); 
  float v_fl                    = fl.getVal();
  float e_fl                    = fl.getError(); //sqrt( (v_nsig1*e_nsig0)*(v_nsig1*e_nsig0)+(v_nsig0*e_nsig1)*(v_nsig0*e_nsig1) )/(v_nsig*v_nsig);
  float v_nqq                   = nqq.getVal();
  float e_nqq                   = nqq.getError();
  float v_nbb                   = nbb.getVal();
  float e_nbb                   = nbb.getError();

  float v_de_qq_cheby           = de_qq_cheby.getVal();
  float e_de_qq_cheby           = de_qq_cheby.getError();

  float v_mbc_qq_argus_par   = mbc_qq_argus_par.getVal();
  float e_mbc_qq_argus_par   = mbc_qq_argus_par.getError();

  float v_k0mass_qq_cheby01     = k0mass_qq_cheby01.getVal();
  float e_k0mass_qq_cheby01     = k0mass_qq_cheby01.getError();
  float v_k0mass_qq_bw_mean     = k0mass_qq_bw_mean.getVal();
  float e_k0mass_qq_bw_mean     = k0mass_qq_bw_mean.getError();
  float v_k0mass_qq_bw_wid      = k0mass_qq_bw_wid.getVal();
  float e_k0mass_qq_bw_wid      = k0mass_qq_bw_wid.getError();
  float v_bw_qq_fraction1       = bw_qq_fraction1.getVal();
  float e_bw_qq_fraction1       = bw_qq_fraction1.getError();

  float v_kpmass_qq_cheby01     = kpmass_qq_cheby01.getVal();
  float e_kpmass_qq_cheby01     = kpmass_qq_cheby01.getError();
  float v_kpmass_qq_bw_mean     = kpmass_qq_bw_mean.getVal();
  float e_kpmass_qq_bw_mean     = kpmass_qq_bw_mean.getError();
  float v_kpmass_qq_bw_wid      = kpmass_qq_bw_wid.getVal();
  float e_kpmass_qq_bw_wid      = kpmass_qq_bw_wid.getError();
  float v_bw_qq_fraction2       = bw_qq_fraction2.getVal();
  float e_bw_qq_fraction2       = bw_qq_fraction2.getError();

  float v_nn_tran2_qq_gau1_mean = nn_tran2_qq_gau1_mean.getVal();
  float e_nn_tran2_qq_gau1_mean = nn_tran2_qq_gau1_mean.getError();
  float v_nn_tran2_qq_gau1_wid  = nn_tran2_qq_gau1_wid.getVal();
  float e_nn_tran2_qq_gau1_wid  = nn_tran2_qq_gau1_wid.getError();

  float v_k0hcos_qq_c0  = k0hcos_qq_c0.getVal();
  float v_k0hcos_qq_c1  = k0hcos_qq_c1.getVal();
  float v_k0hcos_qq_c2  = k0hcos_qq_c2.getVal();
  float v_k0hcos_qq_c3  = k0hcos_qq_c3.getVal();

  float v_kphcos_qq_c0  = kphcos_qq_c0.getVal();
  float v_kphcos_qq_c1  = kphcos_qq_c1.getVal();
  float v_kphcos_qq_c2  = kphcos_qq_c2.getVal();
  float v_kphcos_qq_c3  = kphcos_qq_c3.getVal();
  float v_kphcos_qq_c4  = kphcos_qq_c4.getVal();

  float e_k0hcos_qq_c0  = k0hcos_qq_c0.getError();
  float e_k0hcos_qq_c1  = k0hcos_qq_c1.getError();
  float e_k0hcos_qq_c2  = k0hcos_qq_c2.getError();
  float e_k0hcos_qq_c3  = k0hcos_qq_c3.getError();

  float e_kphcos_qq_c0  = kphcos_qq_c0.getError();
  float e_kphcos_qq_c1  = kphcos_qq_c1.getError();
  float e_kphcos_qq_c2  = kphcos_qq_c2.getError();
  float e_kphcos_qq_c3  = kphcos_qq_c3.getError();
  float e_kphcos_qq_c4  = kphcos_qq_c4.getError();

  float v_de_bb_cheby     = de_bb_cheby0.getVal();
  float v_mbc_bb_par      = mbc_bb_argus_par.getVal();
  float v_k0mass_bb_cheby = k0mass_bb_cheby.getVal();
  float v_kpmass_bb_cheby = kpmass_bb_cheby0.getVal();

  float v_nn_tran2_qq_frac = nn_tran2_qq_frac.getVal();
  float e_nn_tran2_qq_frac = nn_tran2_qq_frac.getError();


  result[0] = v_nsig;
  result[1] = e_nsig;
  result[2] = v_fl;
  result[3] = e_fl;
  result[4] = v_nqq;
  result[5] = e_nqq;
  result[6] = v_nbb;
  result[7] = e_nbb;
  result[8] = v_de_qq_cheby;
  result[9] = e_de_qq_cheby;
  result[10] = v_mbc_qq_argus_par;
  result[11] = e_mbc_qq_argus_par;
  result[12] = v_k0mass_qq_cheby01;
  result[13] = e_k0mass_qq_cheby01;

  result[14] = v_k0mass_qq_bw_mean;
  result[15] = e_k0mass_qq_bw_mean;
  result[16] = v_k0mass_qq_bw_wid;
  result[17] = e_k0mass_qq_bw_wid;

  result[18] = v_bw_qq_fraction1;
  result[19] = e_bw_qq_fraction1;

  result[20] = v_kpmass_qq_cheby01;
  result[21] = e_kpmass_qq_cheby01;
  result[22] = v_kpmass_qq_bw_mean;
  result[23] = e_kpmass_qq_bw_mean;
  result[24] = v_kpmass_qq_bw_wid;
  result[25] = e_kpmass_qq_bw_wid;

  result[26] = v_bw_qq_fraction2;
  result[27] = e_bw_qq_fraction2;
  result[28] = v_nn_tran2_qq_gau1_mean;
  result[29] = e_nn_tran2_qq_gau1_mean;
  result[30] = v_nn_tran2_qq_gau1_wid;
  result[31] = e_nn_tran2_qq_gau1_wid;

  result[32] = v_k0hcos_qq_c0;
  result[33] = v_k0hcos_qq_c1;
  result[34] = v_k0hcos_qq_c2;
  result[35] = v_k0hcos_qq_c3;

  result[36] = v_kphcos_qq_c0;
  result[37] = v_kphcos_qq_c1;
  result[38] = v_kphcos_qq_c2;
  result[39] = v_kphcos_qq_c3;
  result[40] = v_kphcos_qq_c4;

  result[41] = e_k0hcos_qq_c0;
  result[42] = e_k0hcos_qq_c1;
  result[43] = e_k0hcos_qq_c2;
  result[44] = e_k0hcos_qq_c3;

  result[45] = e_kphcos_qq_c0;
  result[46] = e_kphcos_qq_c1;
  result[47] = e_kphcos_qq_c2;
  result[48] = e_kphcos_qq_c3;
  result[49] = e_kphcos_qq_c4;

  result[50] = v_de_bb_cheby;
  result[51] = v_mbc_bb_par;
  result[52] = v_k0mass_bb_cheby;
  result[53] = v_kpmass_bb_cheby;

  result[54] = v_nn_tran2_qq_frac;
  result[55] = e_nn_tran2_qq_frac;



  delete xxx;
  delete yyy01;
  delete yyy02;
  delete yyy03;
  delete yyy04;
  delete data01;
  delete data02;
  delete final_result;

  delete data_siga;
  delete data_qqa;
  delete data_bba;
  delete data_rarea;
  delete data_sigb;
  delete data_qqb;
  delete data_bbb;
  delete data_rareb;

  delete w11;
  delete w12;
  delete w21;
  delete w22;
  delete w31;
  delete w32;
  delete w41;
  delete w42;
  delete w51;
  delete w52;
  delete w53;
  delete w54;
  delete w55;
  delete w56;
  delete w57;
  delete w58;
  delete w59;
  delete w60;
  delete w61;
  delete w62;
  delete w63;
  delete w64;
  delete w67;

  return;
 
}


//
// converting : float > string
//
string ftos(float f){
   stringstream s;
   s << f;
   return s.str();
}

//
// converting : int > string
//
string itos(int f){
   stringstream s;
   s << f;
   return s.str();
}


//
// roudning up
//
float rounding(float a, int p){
  float tmp = pow(10.0, p);
  float b   = floor(a*tmp+0.5)/tmp;
  return b;
}

//        
// get File Size      
//        
long getFSize(string filename)    
{         
  long file_size = 0; 
  FILE *fp = fopen( filename.c_str() , "r" ); 
  if ( fp == NULL )   
    return -1;        
  fseek( fp, 0l , SEEK_END );     
  file_size = ftell(fp);          
  fclose(fp);         
  return file_size;   
}

//
//******************************** For relativistic BW *************************************
//

ClassImp(RooRelativisticBreitWigner) 

RooRelativisticBreitWigner::RooRelativisticBreitWigner(const char *name, const char *title,
		RooAbsReal& _m,
		RooAbsReal& _mass_J,
		RooAbsReal& _width_J,
		RooAbsReal& _r,
		RooAbsReal& _J,
		double _m1,
		double _m2) :
		RooAbsPdf(name,title),
		m("m","m",this,_m),
		mass_J("mass_J","mass_J",this,_mass_J),
		width_J("width_J","width_J",this,_width_J),
		r("r","r",this,_r),
		J("J","J",this,_J),
		m1(_m1),
		m2(_m2) {
}


RooRelativisticBreitWigner::RooRelativisticBreitWigner(const RooRelativisticBreitWigner& other, const char* name) :
	   RooAbsPdf(other,name),
	   m("m",this,other.m),
	   mass_J("mass_J",this,other.mass_J),
	   width_J("width_J",this,other.width_J),
	   r("r",this,other.r),
	   J("J",this,other.J),
	   m1(other.m1),
	   m2(other.m2) {
}


double RooRelativisticBreitWigner::CalculateMomentum(double mm) const {
	const double add_12 = m1 + m2;
	const double sub_12 = m1 - m2;
	return std::sqrt((mm*mm - add_12*add_12)*(mm*mm - sub_12*sub_12))/(2.0*mm);
}


Double_t RooRelativisticBreitWigner::evaluate() const {
	double q = CalculateMomentum(m);
	double q_J = CalculateMomentum(mass_J);
	if(std::isnan(q) || std::isnan(q_J)) return 0.0;
	std::complex<double> amplitude = ComplexPDF::RelativisticBreitWigner(m,mass_J,width_J,r,(int)J,q,q_J);
	return std::norm(amplitude);
}

namespace ComplexPDF {

  double CalculateMomentum(const double m, const double m1, const double m2) {
    double add_12 = m1 + m2;
    double sub_12 = m1 - m2;
    return std::sqrt((m*m - add_12*add_12)*(m*m - sub_12*sub_12))/(2.0*m);
  }

  double BlattWeisskopfFormFactor(const double q, const double q_J, const double r, const unsigned int J) {
    double r2 = r*r;
    double q2 = q*q;
    double q_J2 = q_J*q_J;
    if(J == 0) { return 1.0; } 
    else if (J == 1) { return std::sqrt((1 + r2*q_J2) / (1 + r2*q2)); } 
    else if (J == 2) { return std::sqrt((9 + 3*r2*q_J2 + r2*r2*q_J2*q_J2) / (9 + 3*r2*q2 + r2*r2*q2*q2)); } 
    else {
      std::cerr << "ERROR: Spin " << J << " not implemented for BlattWeisskopfFormFactor(..) in libComplexPDF" << std::endl;
      std::exit(EXIT_FAILURE);
    }
  }

  std::complex<double> RelativisticBreitWigner(const double m, const double mass_J, const double width_J, const double r, const unsigned int J, const double q, const double q_J) {
    double q_ratio = q / q_J;
    double m_ratio = mass_J / m;
    double mass_J2 = mass_J*mass_J;
    double m2 = m*m;
    double delta_J = 0;
    if(J == 0) {
      double width = width_J * ( m_ratio ) * ( q_ratio );
      delta_J = std::atan( ( mass_J * width ) / ( mass_J2 - m2 ) );
    } else if (J == 1) {
      double width = width_J * ( m_ratio ) * ( BlattWeisskopfFormFactor(q,q_J,r,1) * BlattWeisskopfFormFactor(q,q_J,r,1) ) * std::pow( q_ratio, 3 );
      delta_J = std::atan( ( mass_J * width ) / ( mass_J2 - m2 ) );
    } else if(J == 2) {
      double width = width_J * ( m_ratio ) * ( BlattWeisskopfFormFactor(q,q_J,r,2) * BlattWeisskopfFormFactor(q,q_J,r,2) ) * std::pow( q_ratio, 5 );
      delta_J = std::atan( ( mass_J * width ) / ( mass_J2 - m2 ) );
    } else {
      std::cerr << "ERROR: Spin " << J << " not implemented for RelativisticSpinBreitWigner(..) in libComplexPDF" << std::endl;
      std::exit(EXIT_FAILURE);
    }
    double phsp = m / q;
    std::complex<double> amplitude = std::polar(std::sin(delta_J), delta_J);
    return phsp * amplitude;
  }

} // ComplexPDF

//
//******************************************************************************************
//
